Library and functions

############
##Libraries#
############

suppressMessages(library(plot.matrix))      #plot correlation matrix
suppressMessages(library(RColorBrewer))     #color palette
suppressMessages(library(R.matlab))         #import mat files
suppressMessages(library(dplyr))            #mean and sd by group
suppressMessages(library(splitstackshape))  #random stratified sampling
suppressMessages(library(reshape2))         #plotting libraries
suppressMessages(library(ggplot2))
suppressMessages(library(gridExtra))
suppressMessages(library(grid))
suppressMessages(library(egg))
suppressMessages(library(pbapply))          #progress bar in apply functions

Load and preprocess data

#Choose rat and session
rat <- "rat08"; folder <- "190410"; wind <- 128

########################
#Load all tetrode files#
########################

setwd(paste("D:/Sergio/MSI Prestige/PhD/GCM_rats/IPHYS_David/data/", rat, "/", rat, "_", folder,"_statistics", sep=""))
listFrames <- list.files(pattern="spikesFrame.*?\\.csv")
spikesTetr <- lapply(listFrames, read.csv)

#Get number of units detected in every tetrode
tetrUnits <- sapply(spikesTetr, function(x) length(unique(x$unit)))

#Change the unit number in the full data frame
sumIni <- c(0, cumsum(tetrUnits)[-length(tetrUnits)])
for (i in 1:length(tetrUnits)) {
  spikesTetr[[i]]$unit <- spikesTetr[[i]]$unit + sumIni[i]
}

##########################
#Load noise removal times#
##########################
setwd("D:/Sergio/MSI Prestige/PhD/GCM_rats/IPHYS_David/data/noise_removal")
noiseRem <- readMat(paste(rat, "_", folder, "_window_", wind, ".mat", sep=""))[[1]]
#1=noise, 0=no noise
noiseThres <- floor(0.8*length(spikesTetr))
int2rem <- colSums(noiseRem) >= noiseThres

Remove noise

Function for all tetrodes

#Function to label each spike as noise ot not
noiseMark <- function(events, noiseMat, wind) {
  #1=noise, 0=no noise
  int2rem <- colSums(noiseMat) >= noiseThres
  #Time interval vector
  timesInt <- (0:(dim(noiseMat)[2]+1))*wind
  #Find to which interval corresponds each spike event
  whichInt <- sapply(events$spikes*1e3, function(x) findInterval(x, timesInt))
  #Label each spike as noise or not
  noise <- sapply(whichInt, function(x) as.integer(int2rem[x]))
  #Store in data frame
  spikes2denoise <- data.frame(spikes = events$spikes, unit = events$unit, stimuli = events$stimuli, interval = whichInt, noise = noise)
  spikes2denoise <- spikes2denoise[complete.cases(spikes2denoise), ]
  return(spikes2denoise)
}

# Store all tetrodes in a list
list2denoise <- pblapply(spikesTetr, function(x) noiseMark(x, noiseRem, wind))
saveRDS(list2denoise,
        file=paste("D:/Sergio/MSI Prestige/PhD/GCM_rats/IPHYS_David/data/noise_removal/", rat, "_", folder,
                   "_NoiseRemoval_window", wind, ".RData", sep=""))
# list2denoise <- readRDS(paste("D:/Sergio/MSI Prestige/PhD/GCM_rats/IPHYS_David/data/noise_removal/", rat, "_", folder,
#                               "_NoiseRemoval_window", wind, ".RData", sep=""))

#Check percentages
for (i in 1:dim(noiseRem)[1]) {
  print(table(list2denoise[[i]]$noise))
}

    0     1 
 2677 18017 

   0    1 
7153 9687 

    0     1 
 1527 18275 

    0     1 
14948 17581 

    0     1 
 1266 16883 

    0     1 
 7078 19308 

   0    1 
7652 8106 

    0     1 
 8471 17496

Plot clustering without noise

tetrs <- 1:dim(noiseRem)[1] #Rat01, 02, 08
threshold <- c("1", "005", "01", "2", "1", "01", "05", "2") #Rat08 190410
# tetrs <- c(1:5, 8) #Rat 03
# threshold <- c("005", "1", "05", "01", "005", "4") #Rat03 190306
# threshold <- c("2", "01", "005", "015", "03", "005", "05") #Rat02 190306
# threshold <- c("05", "2", "025", "025", "10", "2", "8", "3") #Rat01 190308

You must enable Javascript to view this page properly.

Correlation matrix and raster

Preprocess data

#Join all tetrodes into one data frame 
spikesFull <- dplyr::bind_rows(list2denoise); spikesFull <- spikesFull[order(spikesFull$unit), ]

#Stimuli trains
trains <- readMat(paste("D:/Sergio/MSI Prestige/PhD/GCM_rats/IPHYS_David/data/", rat, "/Trains_", rat, "_", folder, 
                         "_RAW.mat", sep=""))
feeder <- unlist(trains[[1]][1][[1]][2])
lever <- unlist(trains[[1]][1][[1]][3]); leverSec <- rep(c("press","release"), length(lever)/2)
#CAUTION: IF DATA IS "MERGED", SOFTWARE IS THE 5TH ELEMENT, IF NOT, 4TH
software <- unlist(trains[[1]][1][[1]][4])
#0 begin/end baseline, 1 reward, 2 non reward left, 3 non reward right, 4 end behavioral
stimOrder <- c(0,0,1,2,3,2,1,3,2,1,3,1,3,2,3,1,2,3,1,2,1,2,3,2,1,3,2,1,3,
               1,3,2,3,1,2,3,1,2,1,2,3,2,1,3,2,1,3,1,3,2,3,1,2,3,1,2,4)
stimCount <- c(1:57)

Compute correlation matrix

countsFinal <- list()

#Parameters
binSize <- 2
freq <- 1e4
spikesNoNoise <- spikesFull[spikesFull$noise == 0, ]
maxTime <- floor(max(spikesNoNoise$spikes)*freq)

for (g in 1:3) {
  #Subset stimuli
  spikeSubset <- spikesNoNoise[spikesNoNoise$stimuli == g, ]
  minSpike <- min(spikeSubset$spikes); maxSpike <- max(spikeSubset$spikes)
  
  #Generate full time vector with freq precision
  vecTimes <- seq(0, maxTime, by = 1)
  stim2 <- floor(software[stimOrder == g]*freq); stim2end <- stim2 + 30*freq
  
  #Storage objects
  countsToCorr <- list()  #for spike counting

  ####COUNT SPIKES####
  for (i in unique(spikesNoNoise$unit)) {
    #Spikes time to freq precision
    freqSpikes <- floor(spikeSubset[spikeSubset$unit == i, 1]*freq)
    #Spike detection in full time vector
    spikesPresence <- vecTimes %in% freqSpikes
    
    ##
    stim2Count <- vector()
    for (j in 1:length(stim2)) {
      stim2Count <- c(stim2Count, spikesPresence[stim2[j]:stim2end[j]])
    }
    ##
    
    ##
    sumInt <- seq(1, length(stim2Count), binSize*freq)
    binCounts <- vector()
    for (k in sumInt[-length(sumInt)]) {
      binCounts <- c(binCounts, sum(stim2Count[k:(k+(binSize*freq))]))
    }
    ##
    countsToCorr[[i]] <- binCounts
  }
  
  #Store in one matrix
  lenMat <- length(countsToCorr[[1]])
  countsToCorr.matrix <- matrix(unlist(countsToCorr), ncol = lenMat, byrow=T)
  countsToCorr.matrix[is.na(countsToCorr.matrix)] <- 0
  countsFinal[[g]] <- countsToCorr.matrix
  
}

uu <- countsFinal[[1]]; vv <- countsFinal[[2]]; yy <- countsFinal[[3]]; zz <- cbind(uu,vv, yy)

corMat <- cor(zz)
spear.corMat <- cor(zz, method=c("spearman"))
kend.corMat <- cor(zz, method=c("kendall"))

Plot raster plot

####SPIKE COUNTS PLOT######
rastermat <- cbind(uu, vv, yy)
rastermat[rastermat > 30] <- 30 #adjust for every rat and session
melted_rastMat <- melt(rastermat)
midpoint <- (min(rastermat)+max(rastermat))/2
meanpoint <- mean(rastermat)
medianpoint <- median(rastermat)

ggplot(data = melted_rastMat, aes(x=Var2, y=Var1, fill=value)) + 
  ggtitle("Spike count") +
  geom_tile() +
  scale_fill_gradient2(low = "red", high = "yellow", mid = "white", midpoint = midpoint,  
                       limit = c(min(rastermat),max(rastermat)), space = "Lab", name="Spikes count") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(),
        plot.title = element_text(hjust = 0.5, size = (9))) + 
  scale_y_continuous(trans = "reverse") +
  # geom_hline(yintercept=c(5.5,8.5,10.5,12.5,15.5,17.5), linetype="solid", color = "green",  size=1) +  
  geom_vline(xintercept=c(dim(uu)[2]+0.5, 2*dim(uu)[2]+0.5), linetype="solid", color = "blue", size=2)

Repeated measures correlation

colfunc.shape <- colorRampPalette(c("navy", "white", "red"))
numcols <- 100
colcol.shape <- colfunc.shape(numcols)
plot(rep(1,numcols),col=(colcol.shape), pch=19,cex=2)

nunits <- dim(uu)[1]

#Format stim 1 mat
binSize <- 2
numtot <- 18*(30/binSize)
sepVec <- seq(1, numtot + 5, 30/binSize); sepEnd <- sepVec[-1] - 1 
sepList <- list()
for (i in 1:length(sepEnd)) {
  sepList[[i]] <- uu[, sepVec[i]:sepEnd[i]]
}

uu.RealStim <- do.call(rbind, sepList)
uu.Neurons <- rep(1:nunits, 18)
uu.toRM <- data.frame(unit = factor(uu.Neurons), uu.RealStim)

# aa <- rmcorr(unit, "X1", "X2", uu.toRM)

#Format stim 2 mat
sepList <- list()
for (i in 1:length(sepEnd)) {
  sepList[[i]] <- vv[, sepVec[i]:sepEnd[i]]
}

vv.RealStim <- do.call(rbind, sepList)
vv.Neurons <- rep(1:nunits, 18)
vv.toRM <- data.frame(unit = factor(vv.Neurons), vv.RealStim)

#Format stim 3 mat
sepList <- list()
for (i in 1:length(sepEnd)) {
  sepList[[i]] <- yy[, sepVec[i]:sepEnd[i]]
}

yy.RealStim <- do.call(rbind, sepList)
yy.Neurons <- rep(1:nunits, 18)
yy.toRM <- data.frame(unit = factor(yy.Neurons), yy.RealStim)

#################

Mean by neuron

uu.toRMmeans <- matrix(unlist(by(uu.toRM[,-1], uu.toRM$unit, colMeans)), nunits, (30/binSize), byrow = TRUE)
vv.toRMmeans <- matrix(unlist(by(vv.toRM[,-1], vv.toRM$unit, colMeans)), nunits, (30/binSize), byrow = TRUE)
yy.toRMmeans <- matrix(unlist(by(yy.toRM[,-1], yy.toRM$unit, colMeans)), nunits, (30/binSize), byrow = TRUE)

zz.toRMmeans <- cbind(uu.toRMmeans, vv.toRMmeans, yy.toRMmeans)

corMat.toRMmeans <- cor(zz.toRMmeans, method="kendall")

meancolDiagKen <- corMat.toRMmeans
diag(meancolDiagKen) <- 0

Ggplot

library(reshape2)
library(ggplot2)

melted_cormat <- melt(meancolDiagKen)

ggplot(data = melted_cormat, aes(x=Var1, y=Var2, fill=value)) + 
  geom_tile() + 
  scale_fill_gradient2(low = "navy", high = "red", mid = "white", midpoint = (min(meancolDiagKen)+max(meancolDiagKen))/2,  
                       limit = c(min(meancolDiagKen),max(meancolDiagKen)), space = "Lab", name="Kendall\nCorrelation") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, size = 12, hjust = 1)) +
  coord_fixed() +
  geom_hline(yintercept=c((30/binSize)+0.5, (60/binSize)+0.5), linetype="solid", color = "green", size=1) + 
  geom_vline(xintercept=c((30/binSize)+0.5, (60/binSize)+0.5), linetype="solid", color = "green", size=1)

Format data for means

#Comparison labels
division1 <- dim(uu.toRMmeans)[2]
division2 <- dim(uu.toRMmeans)[2] + dim(vv.toRMmeans)[2]
division3 <- dim(uu.toRMmeans)[2] + dim(vv.toRMmeans)[2] + dim(yy.toRMmeans)[2];

labs <- c(rep(c(rep("s1s1", dim(uu.toRMmeans)[2]), rep("s1s2", dim(uu.toRMmeans)[2]), rep("s1s3", dim(yy.toRMmeans)[2])), division1),
          rep(c(rep("s2s1", dim(vv.toRMmeans)[2]), rep("s2s2", dim(vv.toRMmeans)[2]), rep("s2s3", dim(yy.toRMmeans)[2])), division1),
          rep(c(rep("s3s1", dim(vv.toRMmeans)[2]), rep("s2s3", dim(vv.toRMmeans)[2]), rep("s3s3", dim(yy.toRMmeans)[2])), division1))
labMat <- matrix(labs, division3, division3, byrow = FALSE)
halfLabs <- labMat[lower.tri(labMat,  diag=FALSE)]
halfCorMat <- corMat.toRMmeans[lower.tri(corMat.toRMmeans, diag=FALSE)]

#Mean, SE and SE
fullComp.means <- by(halfCorMat, halfLabs, mean)
fullComp.sd <- by(halfCorMat, halfLabs, sd)
fullComp.se <- by(halfCorMat, halfLabs, function(x) sd(x)/sqrt(length(x)))

#CI95
alpha <- 0.05
t <- qt((1-alpha)/2 + .5, 108-1) 
fullComp.CI95 <- t*fullComp.se
#All in data frame
fullComp.frame <- data.frame(Comp = unique(halfLabs), Mean = as.vector(fullComp.means), 
                             SD = as.vector(fullComp.sd), SE = as.vector(fullComp.se), CI95 = as.vector(fullComp.CI95))

#Barplot
library(ggplot2)
ggplot(fullComp.frame) +
  geom_bar( aes(x=Comp, y=Mean), stat="identity", fill="forestgreen", alpha=0.5, width = 0.5) +
  geom_errorbar( aes(x=Comp, ymin=Mean-SE, ymax=Mean+SE), width=0.2, colour="orange", alpha=0.9, size=1) +
  ggtitle("using SE") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(),
          axis.line = element_line(colour = "black"))

Stats

# dataStat <- data.frame(Comp = labs,
#                        Corrs = as.vector(corMat[, 1:division1]))

library(DescTools)
dataStat <- data.frame(Comp = halfLabs, Corrs = FisherZ(halfCorMat))

## parametric
summary(aov(Corrs ~ Comp, dataStat))
             Df Sum Sq Mean Sq F value Pr(>F)    
Comp          5  85.80  17.160   362.2 <2e-16 ***
Residuals   984  46.62   0.047                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pairwise.t.test(dataStat$Corrs, dataStat$Comp, p.adj = "bonf")

    Pairwise comparisons using t tests with pooled SD 

data:  dataStat$Corrs and dataStat$Comp 

     s1s1    s1s2    s1s3    s2s2    s2s3   
s1s2 < 2e-16 -       -       -       -      
s1s3 < 2e-16 1.00    -       -       -      
s2s2 1.00    < 2e-16 < 2e-16 -       -      
s2s3 0.88    < 2e-16 < 2e-16 0.61    -      
s3s3 6.2e-07 < 2e-16 < 2e-16 2.9e-07 8.5e-05

P value adjustment method: bonferroni 
## non parametric
kruskal.test(Corrs ~ Comp, dataStat)

    Kruskal-Wallis rank sum test

data:  Corrs by Comp
Kruskal-Wallis chi-squared = 643.12, df = 5, p-value < 2.2e-16
pairwise.wilcox.test(dataStat$Corrs, dataStat$Comp,  p.adj = "bonf")

    Pairwise comparisons using Wilcoxon rank sum test 

data:  dataStat$Corrs and dataStat$Comp 

     s1s1    s1s2    s1s3    s2s2    s2s3   
s1s2 < 2e-16 -       -       -       -      
s1s3 < 2e-16 1.00000 -       -       -      
s2s2 1.00000 < 2e-16 < 2e-16 -       -      
s2s3 1.00000 < 2e-16 < 2e-16 1.00000 -      
s3s3 0.00034 < 2e-16 < 2e-16 0.00043 0.00177

P value adjustment method: bonferroni 
## stratified random sampling
stratDat <- stratified(dataStat, "Comp", 1000)
stratDat %>% group_by(Comp) %>% summarise(comp_mean = mean(Corrs), comp_sd = sd(Corrs))
# A tibble: 6 x 3
  Comp  comp_mean comp_sd
  <fct>     <dbl>   <dbl>
1 s1s1      0.956   0.257
2 s1s2      0.436   0.179
3 s1s3      0.412   0.168
4 s2s2      0.952   0.269
5 s2s3      1.01    0.249
6 s3s3      1.12    0.215
kruskal.test(Corrs ~ Comp, stratDat)

    Kruskal-Wallis rank sum test

data:  Corrs by Comp
Kruskal-Wallis chi-squared = 643.12, df = 5, p-value < 2.2e-16
pairwise.wilcox.test(stratDat$Corrs, stratDat$Comp,  p.adj = "bonf")

    Pairwise comparisons using Wilcoxon rank sum test 

data:  stratDat$Corrs and stratDat$Comp 

     s1s1    s1s2    s1s3    s2s2    s2s3   
s1s2 < 2e-16 -       -       -       -      
s1s3 < 2e-16 1.00000 -       -       -      
s2s2 1.00000 < 2e-16 < 2e-16 -       -      
s2s3 1.00000 < 2e-16 < 2e-16 1.00000 -      
s3s3 0.00034 < 2e-16 < 2e-16 0.00043 0.00177

P value adjustment method: bonferroni